program example_FEM
	use MySparseOperator
	use printVecMtx
	use lFEMRect

	implicit none
	integer, dimension(:), allocatable :: ifa,jfa,ifb,jfb
	real(kind=8), dimension(:), allocatable :: fa,fb, Fv, uh
	integer :: mtype, m, n
	real(kind=8) :: derrs(4)


! get the stiffness matrix on the Richardarsion mesh of [0,1]x[0,1], all functions are valued by double real
 	if ( allocated(fa) ) deallocate(fa) 
 	if ( allocated(ifa) ) deallocate(ifa)
 	if ( allocated(jfa) ) deallocate(jfa)
 	if ( allocated(fb) ) deallocate(fb)
 	if ( allocated(ifb) ) deallocate(ifb)
 	if ( allocated(jfb) ) deallocate(jfb)
 	call lfemrect_0Dirichlet_dstiffmtx( (/0.d0,0.d0,1.d0,1.d0/),4,5,c_fun,a_fun,fa,ifa,jfa,fb,ifb,jfb)
 	write(*,*) '=============Stiffness Matrix  A ============================='
 	call prtVec(fa,'f8.3')
 	call prtVec(ifa,'I4')
 	call prtVec(jfa,'I4')
 	write(*,*) '=============Stiffness Matrix  B ============================='
	call prtVec(fb,'f8.3')
 	call prtVec(ifb,'I4')
 	call prtVec(jfb,'I4')
 	write(*,*) '=============Stiffness Matrix  A(1:8,1:8) ============================='
 	call mso_prtspmtxr(fa, ifa, jfa, (/1,9/), (/1,9/), 'f10.2')

 	if ( allocated(fa) ) deallocate(fa) 
 	if ( allocated(ifa) ) deallocate(ifa)
 	if ( allocated(jfa) ) deallocate(jfa)
 	call lfemrect_0Dirichlet_dstiffmtx( (/0.d0,0.d0,1.d0,1.d0/),4,3,c_fun,a_fun,fa,ifa,jfa)
 	write(*,*) '=============Stiffness Matrix  A without boundary ============================='
 	call prtVec(fa,'f8.3')
 	call prtVec(ifa,'I4')
 	call prtVec(jfa,'I4')

!  get the right vector on the Richardarsion mesh of [0,1]x[0,1], all functions are valued double real
	if ( allocated(Fv) ) deallocate(Fv) 
	call lfemrect_0Dirichlet_drightvec( (/0.d0,0.d0,1.d0,1.d0/),2,2,f_fun, Fv)
	write(*,*) '=============right Vector F ============================='
	call prtVec(Fv, 'f8.3')
	if ( allocated(Fv) ) deallocate(Fv) 
	call lfemrect_0Dirichlet_drightvec( (/0.d0,0.d0,1.d0,1.d0/),4,3,f_fun, Fv,'n')
	write(*,*) '=============right Vector F without boundary ============================='
	call prtVec(Fv, 'f8.3')

!  get the errors with zeros boundary condition on Cartesian mesh of [0,1]x[0,1]
	mtype = 2	! the stiffness matrix is symmetric definite for this problem
	do m = 3, 18, 3
		n = m
		! get the linear fem solution with zeros boundary condition on Cartesian mesh of [0,1]x[0,1]
		if ( allocated(uh) ) deallocate(uh)
		call lfemrect_0Dirichlet_dsolution( (/0.d0,0.d0,1.d0,1.d0/),m,n,c_fun,a_fun,f_fun,uh,mtype)
		call lfemrect_0Dirichlet_derrors( (/1,1,1,1/),(/0.d0,0.d0,1.d0,1.d0/),m,n,u_fun,uh,derrs)
		write(*,'(I10,I10,4d15.3)') m,n,derrs
	end do

 	if ( allocated(fa) ) deallocate(fa) 
 	if ( allocated(ifa) ) deallocate(ifa)
 	if ( allocated(jfa) ) deallocate(jfa)
 	if ( allocated(fb) ) deallocate(fb)
 	if ( allocated(ifb) ) deallocate(ifb)
 	if ( allocated(jfb) ) deallocate(jfb)
 	if ( allocated(Fv) ) deallocate(Fv)
 	if ( allocated(uh) ) deallocate(uh) 
 contains
 	subroutine c_fun(x,y,v1,v2)
 		real(kind=8), dimension(:), intent(in) :: x,y
 		real(kind=8), dimension(:), intent(out) :: v1,v2
 		v1 = 1
 		v2 = 1
 	end subroutine c_fun
 	subroutine a_fun(x,y,v)
 		real(kind=8), dimension(:), intent(in) :: x,y
 		real(kind=8), dimension(:), intent(out) :: v
 		v = 1
 	end subroutine a_fun
 	subroutine f_fun(x,y,v)
 		real(kind=8), dimension(:), intent(in) :: x,y
 		real(kind=8), dimension(:), intent(out) :: v
 		v = 2*( x*(1-x) + y*(1-y) ) + x*(1-x)*y*(1-y)
 	end subroutine f_fun
 	subroutine u_fun(x,y,v,vx,vy)
 		real(kind=8), dimension(:,:), intent(in) :: x,y
		real(kind=8), dimension(:,:), intent(out) :: v, vx, vy
		if ( .not.all( shape(x).eq.shape(y) ) ) then
			write(*,*) 'the shape of x does not agree with y in u_fun'
			stop 1
		end if 
		if ( .not.all( shape(x).eq.shape(v) ) ) then
			v = 0.d0
		else
			v = x*(1.d0-x)*y*(1.d0-y)
		end if
		if ( .not.all( shape(x).eq.shape(vx) ) ) then
			vx = 0.d0
		else
			vx = (1.d0-2.d0*x)*y*(1.d0-y)
		end if
		if ( .not.all( shape(x).eq.shape(vy) ) ) then
			vy = 0.d0
		else
			vy = x*(1.d0-x)*(1.d0-2.d0*y)
		end if
	end subroutine u_fun

end program example_FEM